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We present a review of recent developments of simulations of the Vlasov-Maxwell system of 
equations using a Fourier transform method in velocity space. In this method, the distribution 
functions for electrons and ions are Fourier transformed in velocity space, and the resulting set 
of equations are solved numerically. In the original Vlasov equation, phase mixing may lead to 
an oscillatory behavior and sharp gradients of the distribution function in velocity space, which is 
problematic in simulations where it can lead to unphysical electric fields and instabilities and to the 
recurrence effect where parts of the initial condition recur in the simulation. The particle distribution 
function is in general smoother in the Fourier transformed velocity space, which is desirable for the 
numerical approximations. By designing outflow boundary conditions in the Fourier transformed 
velocity space, the highest oscillating terms are allowed to propagate out through the boundary and 
are removed from the calculations, thereby strongly reducing the numerical recurrence effect. The 
outflow boundary conditions in higher dimensions including electromagnetic effects are discussed. 
The Fourier transform method is also suitable to solve the Fourier transformed Wigner equation, 
which is the quantum mechanical analogue of the Vlasov equation for classical particles. 
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I. INTRODUCTION 

The Vlasov equation governs the dynamics of the distribution function of charged particles (electrons, ions) in the 
six-dimensional phase space, consisting of 3 velocity (or momentum) dimensions and 3 position dimensions, plus time. 
It offers an accurate description of a plasma in the coUisionlcss limit, i.e., when the particles are affected by long-range 
electric and magnetic fields only, and when short-range fields from its nearest neighbors can be neglected. 

The most common method to solve the Vlasov equation numerically is the Particle-In-Cell (PIC) method (Birdsall 
and Langdon, 1991; Matsumoto and Omura, 1993), where the Vlasov equation is solved by following the trajectories of 
a set of statistically distributed super-particles, which resolves the particle distribution functions in phase space. Each 
super-particle represents a large number of real particles. PIC simulations have proven to be extremely successful due 
to their relative simplicity and adaptivity. However, the statistical noise of PIC simulations sometimes overshadows 
the physical results, and for some problems, the low-density velocity tail of the particle distribution cannot be resolved 
with high enough accuracy by the super-particles. 

Grid-based Vlasov solvers treat the particle distribution function as a phase fluid that is represented on a grid in 
both position and velocity (or momentum) space. The advantage with grid-based Vlasov solvers is that they do not 
give rise to statistical noise in the simulations, and that the dynamical range is larger than for PIC methods, so that 
the low-density velocity tail of the particle distribution can be resolved much more accurately. A disadvantage with 
grid-based Vlasov solvers in higher dimensions is that the full phase-space has to be represented on a numerical grid 
grid, which makes both the storage of the data in the computer's memory, and the numerical calculations, extremely 
demanding. Another problem is the tendency of the distribution function to become oscillatory in velocity space due 
to phase mixing. This is the cause of Landau damping and other kinetic effects, but may lead to unphysical noise 
and recurrence effects in the numerical solution. This was recognized in early Vlasov simulations and special methods 
were devised to resolve this problem, including the Fourier-Fourier and Hermite expansion methods (Armstrong et 
al, 1970). Various other methods have also been developed for solving the Vlasov equation, the classical and widely 
used time-spitting method (Cheng, 1977; Cheng and Knorr, 1976) where a smoothing operator was applied to remove 
the highest oscillations in velocity space, a Van Leer dissipative scheme (Mangeney et al, 2002), the finite volume 
method (Elkina and Biichner, 2006), a back-substitution method for magnetized plasma (Schmitz and Graucr, 2006), 
etc. Several Eulerian grid-based solvers are reviewed and compared by Arber and Vann (2002) and Filbet and 
Sonnendriicker (2003). 

One special category of methods are transform methods (Armstrong et al., 1970), where the particle distribution 
function is expressed as a sum or integral of basis functions in velocity space. The structuring of velocity space in this 
case leads to the excitation of higher and higher modes in the transformed velocity space, and special care must be taken 
when these excitations reach the highest mode represented in the numerical simulation. For example, for methods 
using Hermite polynomials to resolve the velocity space, the highest-order Hermite polynomials is designed to absorb 
the oscillations in velocity space (Gibclli and Shizgal, 2006; Knorr and Shoucri, 1974), thereby reducing numerical 
recurrence effects strongly. Klimas (1987) and Klimas (1994) devised filtering methods for the Fourier-Fourier method 
to remove the highest oscillations in velocity space. (In the Fourier-Fourier method, the Vlasov equation is Fourier 
transformed both in velocity space and position space, and the resulting equation is solved numerically.) For the 
Fourier method in one, two and three dimensions, Eliasson (2001, 2002, 2003, 2007) designed absorbing boundary 
conditions at the largest Fourier mode in velocity space so that the highest oscillations in velocity space were removed 
from the solution. In this method, the Vlasov equation is Fourier transformed only in velocity space and not in position 
space, so that the electromagnetic fields enter by multiplications instead of convolutions in the transformed Vlasov 
equation. The Fourier transformed Vlasov equation is also interesting to study in its own respect. Some mathematical 
aspects of the Fourier transformed Vlasov equation are given by Neunzert (1971, 1984), and in their analytic and and 
simulation studies, Sedlacec and Nocera (1992, 2002) presented interpretations of the Landau damping, the time echo 
phenomenon, etc., in terms of imperfectly trapped (leaking) waves in the Fourier transformed velocity space. 

The main topic of this paper is a the properties of the Vlasov equation and the Fourier transformed Vlasov equation. 
The full Vlasov-Maxwell system, and the Vlasov-Maxwell system Fourier transformed in velocity space are presented in 
section II. In section HI, we discuss the properties of the Vlasov equation that makes simulations a challenging task. We 
note that oscillations of the particle distribution function in velocity space corresponds to wave packets in the Fourier 



transformed velocity space. Special attention is put on the absorbing artificial boundary conditions in the Fourier 
transformed velocity space, where the highest Fourier modes are absorbed and removed from the calculations. These 
boundary conditions have the following attractive features: (i) They substantially reduce the unphysical reflections at 
the artificial boundaries, thereby reducing unphysical noise and recurrence effects in simulations, (ii) The boundary 
conditions are local in time and involve only boundary points (but they are non-local along the boundary), (iii) The 
boundary conditions together with the interior differential equation defines a well-posed problem. Some examples of 
simulations of one-dimensional kinetic structures, electron and ion holes, are also discussed in section III. In section 
IV, we will briefly mention the discrete approximations used in the numerical simulations of the Fourier transformed 
Vlasov equation. Especially the important topics of the representation of particle velocities by the numerical scheme 
and how to choose domain and grid sizes in the Fourier transformed velocity space are discussed. In sections V and VI, 
we present the generalizations of the Fourier technique to two and three dimensions, respectively. Here the treatment 
of the spatially varying magnetic fleld has to be treated carefully in the design of absorbing boundary conditions int 
the Fourier transformed velocity space. We mention here that the well-posedness of the boundary conditions for the 
one-, two- and three-dimensional Fourier transformed Vlasov equation have been proved by energy estimates (Eliasson, 
2001, 2002, 2007), where an energy norm is non-increasing in time. Finally, in chapter VII, we discuss extensions of 
the Vlasov equation to incorporate quantum and relativistic effects. The quantum analogue to the Vlasov equation 
is the Wigner equation, and we will see that the absorbing boundary conditions used for the Vlasov equation can 
be applied unchanged to the Wigner equation. For the relativistic Vlasov equation, the relativistic gamma factor 
enters into the Vlasov equation and leads to a convolution in the Fourier transformed velocity space. This convolution 
operator is non-local in space and may lead to non-local absorbing boundary conditions in space and time. 

II. THE VLASOV-MAXWELL SYSTEM OF EQUATIONS 

The non-relativistic Vlasov equation 

^+v.V,/„ + ^.V„/„ = (1) 

where the Lorentz force is 

F„ = g„[E + vx (B + Bext)] (2) 

describes the evolution of the distribution function /„ of electrically charged particles of type a (e.g., "electrons" or 
"singly ionized oxygen ions"), each particle having the electric charge qa and mass nia- Here, the magnetic fleld is 
separated into two parts, where Boxt is an external magnetic fleld (e.g., the Earth's geomagnetic fleld), and B is the 
self-consistent part of the magnetic fleld, created by the plasma. One Vlasov equation is needed for each species of 
particles. 

The particles interact via the electromagnetic fleld. The charge and current densities, p and j, act as sources of 
self-consistent electromagnetic flelds according to the Maxwell equations 

V • E = ^ (3) 

V • B = (4) 
VxE.-f (5) 

VxB.,oJ + ^^ (6) 

The charge and current densities are related to the particle number densities n^ and mean velocities v^ as 

a 

and 
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qaUaVa (8) 



respectively, where the particle number densities and mean velocities are obtained as moments of the particle distri- 
bution functions, as 



n„(x,t)= / /„(x,v,t)dS (9) 

and 

1 f°° 

V, (x, t) = -— / v/„ (x, V, t) d^i; (10) 

respectively. 

The Vlasov equation (1) together with the Maxwell equations (3)-(6) and the constitutive equations (7)-(10) form 
a closed, coupled system of nonlinear partial differential equations and integral equations. 

A. The Fourier transformed Vlasov equation 

By using the Fourier transform pair 

/oo 
/„(x,77,t)e-"'-d3r; (11) 

-OO 
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/a(x, r,, t) = j^ / /„(x, V, t)e^^-' d?T^ (12) 



— oo 



the velocity variable v is transformed into a new variable r] and the distribution function /(x,v,t) is changed to a 
new, complex valued, function f{x.,r],t), which obeys the transformed Vlasov equation 

^ - iV, • Vja, - ^(iE • riU + V^ • {[(B + Bext) x r7]/„}) = (13) 

ot nia 

The nabla operators Va; and V^ denote differentiation with respect to x and rj, respectively. 

Equation (13) together with the Maxwell equations (3)-(6) and the constitutive equations (7)-(8) where the particle 
number densities and mean velocities are obtained as 

n^{x,t)^{2TrfU{x,0,t) (14) 

and 

v„(x,i) = -i^^^ [v^7„(x,J7,t)l (15) 

respectively, form a new closed set of equations. One can note that the integrals over infinite v space in Eqs. (9) and 
(10) have been converted to evaluations in rj space in Eqs. (14 and 15). The factor (27r)'^ in Eqs. (12), (14) and (15) 
is valid for three velocity dimensions. For n velocity dimensions this factor is (27r)". 

B. Invariants of the Vlasov-Maxwell system 

The Vlasov equation (1) coupled with (3)-(6) conserves the energy norm 

\\U\?= j j flA\d?x, (16) 

the total number of particles 

Nc.= f f U d^wd^a;, (17) 



the total linear momentuni 



P = /[/vK/.+m./.)d3. + eoExB 



d^x, 



and the total energy 



W = J I ^^^{m,f, + m J,) d^v+UeoE^ 
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The corresponding invariants for the Fourier-transformed Vlasov-Maxwell system (13) and (3)-(6) are: 
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(27r)3 



lU'd'r^d'x, 



N^^ /(27r)3(/„)^^od3 



p = / [-i(27r)3v^(TOi/i + me/c)r,=o + ^oE x B] d^x, 



and 



W = 



- -(2^)3V^(mi/i + me/c)^=o + ^^ EoE' 



B2 

Mo 



d-^x, 



(18) 



(19) 

(20) 
(21) 
(22) 

(23) 



respectively, where the norm (20) follows from (16) via the Parseval relation. These invariants can be used to check 
the accuracy of the numerical scheme. When the system is restricted to a bounded domain in t] space with appropriate 
boundary conditions (discussed in Section VI B below), the norm IL/ajP will be a non- increasing, positive function of 
time, while the other three quantities will still be conserved. 



C. The electrodynamic scalar and vector potentials 

The Maxwell equations (3)"(6) can be written in terms of the scalar and vector potentials <I> and A, which are 
related to the electromagnetic field as 



E= -V$- 



dA 
~dt 

B = Bo + V X A, 



(24) 
(25) 



where Bq is the external magnetic field. Introducing these expressions into Eqs. (3)-(6), and choosing the Lorentz 
condition, 



yields the electrodynamic waves equations 



V-A 



1 ^2$ 

1 a^A 



1 9$ 



V2$ 



= 



£a 



- V^A = Moj. 



(26) 

(27) 
(28) 



In this description, the divergence of the magnetic field is zero, since the divergence of the right hand side of Eq. (25) 
is zero by the vector relation V • (V x A) = 0. The divergence of the electric field can be set to the correct value by 
using the Maxwell equation for the divergence of the electric field (3) together with Eq. (24) , yielding 



V-E = -V2$- V- 



9A 



£o 



(29) 



_V2$^^+V-?^ (30) 

eq at 

This equation for $ conserves the divergence of the electric field. 

By introducing a separate variable F for the time derivative of the vector potential A, the wave equations (28) and 
(30) can be rewritten in a first-order system with respect to time, 

^=r f =c2(V^A + Moj) (31) 

_V2$= — + V-r (32) 

and the electric and magnetic fields are calculated as 

E = -V$-r (33) 

and 

B = V X A (34) 

respectively, in the new variables. 

The system (31)-(34) produces physical electric and magnetic fields regardless of the initial conditions on A and 
r, in the sense that the first two Maxwell equations for the divergences (3)-(4) are fulfilled. Therefore, a consistent 
numerical scheme will also produce physical solutions, up to the local truncation error of the numerical scheme, even 
after a long time; no artificial electric and magnetic charges are created and accumulated by the numerical scheme, 
which could be the case if the two last Maxwell equations (5)-(6) are integrated numerically in time. This general 
property of the system is an advantage, since it is not necessary to use special, divergence-conserving schemes Wagner 
and Schneider (1998) to solve these equations, and it therefore opens up the possibility to switch between different 
numerical methods without having to pay too much attention to the divergences of the electromagnetic field. In 
complicated geometries, it may be a disadvantage that an elliptic equation (32) has to be solved numerically to 
obtain the potential $, while in the simple geometries considered here, with periodic boundary conditions, Eq. (32) 
is efficiently solved by means of Fourier transform techniques. 

D. The reduction of spatial and velocity dimensions 

In the study of problems with certain symmetries, it is sometimes possible to make a choice of the coordinate system 
so that the the problem can be analyzed in a smaller number of dimensions. Numerically, this is very convenient 
because unnecessary information is removed from the problem and a smaller number of sampling points is needed to 
represent the solution on a numerical grid. 

One such assumption is that the problem is homogeneous in one or more dimensions, in which case the derivatives 
in these dimensions vanish. In the study of plane waves in plasma, the number of dimensions in x = (xi,X2, Xa) 
space can be reduced to one dimension, x = (xi,0, 0), so that only derivatives with respect to xi (and not X2 and 
X3) remain. In this manner, the Vlasov equation can be reduced from three spatial and velocity dimensions, to one 
spatial and three velocity dimensions, plus time. 

For the non-relativistic Vlasov equation, it turns out that it is also possible to reduce the number of velocity dimen- 
sions, but in a different manner than for the spatial dimensions. For electrostatic electron waves in an unmagnetized 
plasma, the reduction to one spatial dimension xi also leads to that terms containing factors of V2 and U3, and 
derivatives with respect to V2 and W3, also vanish, giving rise to the system 

^+„,^„^i^^0 (35) 

dt dxi rric dvi 



dEi _ e 

dxi Eq 



— 00 -' —00 ..' —00 



no- / / f{xi,vi,V2,V3,t)dvidv2dv3 



(36) 



The dependence on {v2,V3) only appears in the integral over all velocity space for calculating the electric field Ei. 
Similarly, for waves in a magnetized plasma, propagating in the (a;i,X2) plane perpendicular to the magnetic field 



directed in the X3 direction and with the electric field directed in the (a;i,a:2) plane perpendicular to the magnetic 
field, any dependencies on the distribution in V3 vanish in the Vlasov equation, and one has one Vlasov equation in 
{xi,X2,Vi,V2,t) space for each value on v^. In these cases, it is possible (and convenient) to reduce the number of 
dimensions also in v space. 

In the study of collective phenomena in plasma, the electromagnetic fields do not depend explicitly on the exact 
velocity distribution of particles but on the charge and current densities in x space, calculated as integrals (moments) of 
the distribution function. This makes it possible to reduce the number of velocity dimensions in the Vlasov equation. 
For the case of electrostatic waves in an unmagnetized plasma discussed above, it is simple to show that linear 
combinations of distribution functions with different {v2,V3) are solutions to the one-dimensional Vlasov equation 
(35), because these distribution functions separately are solutions to the same Vlasov equation. In particular, taking 
the limit of a continuous "linear combination," the one-dimensional distribution function 

/OO /"OO 

/ f{xi,vi,V2,V3,t)dv2dv3 (37) 

-00 •/ — 00 

is a solution to the onc-dinicnsional Vlasov equation (35) because the function /(xi,i?i,i'2,i'3, t) is a solution to the 
the Vlasov equation (35) for each value on {v2^v^). The electric field is calculated from Eq. (36) where, by Eq (37), 

/•OO 

fdv,dv2dv3^ / /i°di;i (38) 



— OO J —OO J —OO 



and the resulting one-dimensional Vlasov equation coupled with Gauss' law is 



ot oxi rric ovi 

dEi e 



dxi Eq 



OO 



no - / / (xi,wi,t)dwi 



(40) 



for the unknown function /^^. 

For waves propagating perpendicularly to a magnetic field, mentioned above, it is possible to derive the two- 
dimensional Vlasov-Maxwell system, which depends on the distribution functions in the form 



/ {xi,X2,Vi,V2,t) :^ f{xi,X2,Vi,V2,V3,t)dV3 (41) 

^— OO 

For waves propagating with some angle to the magnetic field, it is more difficult to reduce the number of velocity 
dimensions in the manner described above, since all three velocity components will appear explicitly in the resulting 
Vlasov equation. Even so, the reduction of the number of velocity dimensions for the Vlasov equation has been 
done also for this case, in which a full Vlasov kinetic description is maintained only along one "dominant" spatial 
coordinate, and with the perpendicular dimensions modeled by reduced moment-based methods (Newmann et al., 
2004). 

III. PROPERTIES OF THE VLASOV EQUATION 

A well-known property of the Vlasov equation is that an initially smooth distribution function which evolves in time 
may become increasingly oscillatory in velocity space due to phase mixing of the distribution function. This leads to 
Landau damping and other kinetic effects, but it also makes the numerical solution of the Vlasov equation a challenging 
task (Armstrong et al, 1970; Cheng and Knorr, 1976). We here discuss the main features of the Vlasov equation 
and the numerical difficulties arising from the phase mixing effects. The phase mixing and oscillatory behavior of the 
Vlasov equation in velocity space leads to wave packets in the Fourier transformed velocity space, which points to the 
idea to absorb the highest Fourier modes via boundary conditions in the Fourier transformed Vlasov equation. 



A. The filamentation in velocity space and the numerical recurrence effect 

The behavior of the Vlasov equation is iUustrated by the one-dimensional Vlasov equation for electrons with 
stationary, singly charged ions, coupled with Gauss' law for the electrostatic field. 



dt 
dE _ e 

dx £o 



dx rUe dv 



no 



feix,v,t)dv 



(42) 
(43) 



describing the evolution of the electron distribution function /e in a self-consistent electric field E, where hq is the 
equilibrium electron (and ion) number density, e is the magnitude of the electron charge, me is the electron mass, 
and £o is the electric vacuum permittivity. It can be cast into the dimensionless form 



dt 
dE 

dx 



dfe j^dfe 

V- E^^ = 

dx 



= 1- 



dv 

fe{x,V,t)dv, 



(44) 
(45) 



were the distribution function /e is normalized by nQ/vth,e, time t by wZ.^, space x by r^e, velocity v by vth,e, and the 
electric field E by v^f^ ^rrie/erjje- Here ujpe — (noe^/£o"^e)^'^ is the electron plasma frequency, r^e = (fcsT'e£o/"-oe^)^' ^ 
is the electron Debye length, Wj/j g — [kBTf./mf,)^/'^ is the electron thermal speed, T^ is the electron temperature, and 
ks is Boltzmann's constant. 
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FIG. 1; a) Exact and b) numerical approximation of the electron number density Ue- One can see a numerical recurrence effect 
with periodicity Trccurrcncc = 'i-i^ l(kx Af) ~ 62.8. 



To illustrate the main numerical difficulty of the Vlasov equation, we consider the reduced, free-streaming problem 

,2 



_+,_=0, /(.,«,t = 0) 



—= cos{kxx) exp( — — ) 



which has the exact solution 



/ = {2T:)-^^^Acos[k^{x - vt)] exp(-wV2). 



(46) 



(47) 



This solution becomes more and more oscillatory in the velocity space with increasing time due to the kvt term inside 
the cosinus function. We note that the distribution function / does not decay in time, but, due to phase mixing 



between positive and negative values of /, the number density 

fdv^Aexp{~klt^/2)cos{k^x), (48) 

decays super-exponentially fast with time. Assume now that we have the exact solution / of the electron distribution 
function, and want to calculate the electron number density numerically via a sum representation of the integral over 
velocity space. If we resolve velocity space with an equidistant grid as v = Vj — jAv, j — 0, ±1, ±2, ..., ±M, where 
M is a large integer, then a numerical approximation of the electron number density is 



M 



^ {2'Ky^/'^Acoa[k^{x-jAvt)]exp{-fAv'^/2)Av, 



(49) 



-M 



which turns out to be periodic in time with periodicity Tj-ocurrcncc = 2TT/{k^Av). For example, for kx = 0.5 and 
Av = 0.2, we have Tiecurrence = 27r/(fca; Av) ~ 62.8. While the exact number density decays super-exponentially, we see 
in Fig 1 that in the numerical approximation, the initial condition recurs periodically with periodicity Trccurrcncc ~ 62.8. 
This is the recurrence effect. It is in fact impossible to represent the solution on the grid after a finite time due to 
the Nyquist-Shannon sampling theorem, which states that one needs more than two grid points per wavelength to 
represent the solution an equidistant grid. Here we see from Eq. (47) that the "wavelength" in velocity space is 
At, = 2TT/(kxt). Hence, the sampling theorem Xy/Av > 2 for representing the distribution function on the grid is 
violated for times t > 7r/(fcj,Aw) = Ti-ocurronco/S. 




(a)Thc distribution function f(x, v, t) 
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(b)The real part of the Fourier transformed distribution function 
I(x,r,,t) 



FIG. 2: Phase space plots of the electron distribution function f(x, v, t) and the real part of the Fourier transformed distribution 
function f{x,ri,t) at different times. After Eliasson (2006). 



We now return to the system (44)-(45). By employing the Fourier transform pair 

/oo 
/(a;,77,t)e-""^d77 
-CX3 

the system (44)-(45) is transformed into a new set of equations 



df . d^f 

dt dxdrj 
dE{x,t) 
dx 



iErjf = 



1 -2Trf{x,T],t)rj=o 



(50) 
(51) 

(52) 
(53) 



10 




0.6 40 



30 



t = 70 



W) 




0.20 



O.IB 



0.10 



0.00 



2ir An 



2jr in 



FIG. 3: A closeup of the distribution function f{x,v,t) and of the real part of the Fourier transformed distribution function 
f{x,ri,t) at the time t = 70. After Eliasson (2006). 

for the Fourier transformed distribution function f{x,ri,t). Here / is normalized by no and the Fourier transformed 
velocity variable rj is normalized by v^^^, while the other variables are normalized in the same manner as in Eqs. 
(44)-(45). As initial conditions for Eq. (44), we will use (Armstrong et ai, 1970; Cheng and Knorr, 1976) 



f{x,v,0) = [l + Acos{Kx)]fo{v) 



(54) 



with A = 0.5, kx = 0.5, and foiv) = (27r) ^'^ exp (— u^/2); see Fig. 2(a) at i = 0. The corresponding initial condition 
for the Fourier transformed Vlasov equation (52), shown at i = in Fig. 2(b), is 



fix,r],0) = [l + Acos{k.xx)]fQ{ri), 



(55) 



with fo{ri) — {2tt)^^ exp(— 77^/2). We used the method of Eliasson (2001) to solve the system (44)~(45). The simulation 
domain was set to < x < 47r and < 77 < 120 with N^ ~ 200, N^^ = 600, and the time domain was < t < 70 with 
Nt — 35000 and At ~ 0.002. The numerical dissipation parameter in x space was set to 5 — 0.001. 

As can be seen in Fig. 2(a), the solution has at t = 7 formed filaments with large gradients in velocity space. 
At t = 70, the gradients have further steepened, and two Bernstein-Green-Kruskal (BGK) waves have been formed. 
A closeup of the solution is shown in the left-hand panel of Fig. 3. The initially smooth solution has evolved into 
an oscillatory solution with steep gradients, primarily in v space but also in x space. As a contrast, the Fourier 
transformed solution f{x,ri,t), displayed in Fig. 2(b) for the same times as in Fig. 2(a), does not become oscillatory 
in velocity space. Instead, wave packets are formed and are propagating away from the origin 77 = in the Fourier 
transformed velocity space. 

The smooth solution in the Fourier transformed velocity space is a desirable property for numerical approximations 
of the Vlasov equation. It is possible to give an upper bound on the derivatives in 77 space and on numerical truncation 
errors for numerical schemes that approximate the rj derivatives, which is not possible in the original velocity v space. 
By inspection of the three panels in Fig. 2(a), one can see that the solution has significantly non-zero values only in 
the velocity interval w = — 5 to w = 5. This suggests that the Vlasov equation in the Fourier transformed space has a 
smooth solution. If one assumes that the solution in Fig. 2(a) for all times vanishes as a Gaussian function for large 
V, with the estimate 



\f{x,v,t)\ <Cexp(-7v2) 



(56) 
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for some positive constants C and 7, then the 77 derivatives of the Fourier transformed solution arc bounded as 



dij' 



;fix,il,t) 



= [Use Eq. (51)] = 



< [Use the triangle inequality] < 



1 

2^ 



(if)"e"'V(a;,f,Odw 
\iivre''^-fix,v,t)\dv 



< 



1 

2^ 

1 

2tt 



(57) 



\vr\f{x,v,t)\dv<[ByEq. (56)] 



|w|"Cexp(-7w^)dw = 



1 C 

2^^(«+l)/2' 



where the constant 



[(n-l)/2]!, 



n even 
n odd 



(58) 



and where the symbols ! for the factorial and !! for the semi- factorial have their usual meaning. Thus, by the 
assumption (56) for f{x,v,t) it follows that f{x,ri,t) is infinitely differentiable with respect to rj with the estimate 
(57) for the derivatives. It is therefore possible to make an error estimate of the truncation error of a difference scheme 
used to approximate the 77 derivative in Eq. (52). The 4th-order compact Fade difference scheme, which is used here 
to perform numerical approximations of the first derivatives in r/ space (See Section IV below) has a truncation error 
of size 



|e| < Y^Ary^'max 



a^/ 



drj^ 



where the fifth derivative gives n = 5 in formula (57). This gives the estimate 



e < 



1 A774 C 
2^150^73 



(59) 



(60) 



for the truncation error. It is thus possible to make an error estimate for the numerical differentiation in the Fourier 
transformed velocity 77 space in Eq. (52), which is not possible for a numerical differentiation in the original velocity 
V space in Eq. (44). In the closeup of the solution at time t — 70, displayed in right panel of Fig. 3, a wave packet 
can clearly be seen at ry « 35, which corresponds to the "frequency" of the oscillations in velocity space, seen in the 
left panel of Fig. 3. Since the wave packet is decoupled from the origin 77 = where the electric field is calculated, it 
can be removed from the calculation without immediately affecting the value of the electric field. The wave packets 
eventually reach the artificial boundary ai rj = r]miLx = 120; see the right panel of Fig. 2(b), where they are absorbed, 
as described in the next section. 



B. Outflow boundary conditions in Fourier transformed velocity space 



The idea developed in (Eliasson, 2001, 2002, 2007) is to solve the Vlasov equation in the Fourier transformed velocity 
space and to design absorbing boundary conditions at the largest component; wave packets that reach the artificial 
boundary in the Fourier transformed velocity space are allowed to travel over the boundary and to be removed from 
the solution, while incoming waves are set to zero at the boundary. In this manner a weakly dissipative term is 
introduced in the Vlasov equations, which removes only the highest oscillations in velocity space. By setting the 
artificial boundary further away from the origin in the Fourier transformed velocity space, finer structures is resolved 
in velocity space. 

The idea can be illustrated with the reduced problem 



dl 
dt 



dxdrj 



0. 



A Fourier transform in space (d/dx — > ik) gives the hyperbolic equation 



5/ 
dt 



+ t 



dr] 



= 0, 



(61) 



(62) 
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FIG. 4: The electric field for small-amplitude electrostatic waves, a) and b): outflow boundary conditions in rj space, c): / set 
to zero at boundary rj = rjmax- After Eliasson (2001). 

which has solutions of the form / = foiv ~ k^t). Well-posed outflow boundary conditions for Eq. (62) in 77 space 
are found by setting / to zero at 77 = rimax if ^k < and to zero at 77 = —rjmax if k^ > 0. It can be expressed 
as / = H{kx)f at 77 = rjmax and / = H{—kx)f at 77 = —Tjmax where H is the Heaviside function, here defined as 
H{kx) = for fc^: < and H(kx) = 1 for k^ > 0. Inverse Fourier transforming the boundary conditions, we obtain 
well-posed outflow boundary conditions for /, 



and 



/ = F-i[i?(fc,)F/] at77 = 7;„ 



/ = F-i[iJ(-fc,)F/] at77 = -77„ 



(63) 



(64) 



where F and F~^ are the forward and inverse spatial Fourier transforms. 

It turns out that the outflow boundary condition (63) also works well for the complete Fourier transformed system 
(52)-(53). Figure 4 shows a simulation with a small-amplitude wave in the initial condition so that we have an almost 
linear problem. The wavenumber is set to kx = 0.5 so that the wave is strongly Landau damped. We are comparing 
cases where we have used the outflow boundary condition (63) at 77 = "qmax = 30 [panels a) and b) in Fig. 4] with a 
case where we simply set / to zero at the boundary [panel c)]. [A denser grid in 77 space is used in a) compared to b).] 
In the simulation in Fig. 4, the amplitude of the wave is initially decreasing exponentially as expected from linear 
Vlasov theory. At i « 60ajZ.^, there is a weak recurrence of the waves, and at t w 120w~^ a stronger recurrence takes 
place. We see that the recurrence phenomenon is much weaker for case a) and b) where outflow boundary conditions 
were used in 77 space, while in case c) , where the distribution function was set to zero at the 77 boundary, the amplitude 
of the recurring wave at i = 120 aj~^ is of the same order as in the initial condition. For the case c) the numerical 
results would be useless for a nonlinear problem after the recurrence has taken place at i = 120 wZ.^, while for the 
other cases the linear Landau damping effects are reasonably well represented and a nonlinear problem could be run 
beyond the recurrence time. 




FIG. 5: Large amplitude case. Time development of the energy integral ||/||2 in Eq. (66), normalized by its initial value. 
An interesting question is what is flowing out at the outflow boundary in rj space. It is easy to show that the energy 
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(entropy) integrals ||/|P and ||/|p are conserved in time, where 

ll/ll'-y J fdvdx and WfW^ = J J Ifl'dvdx. (65) 

if periodic boundary conditions are used in space. They are related as ||/|p — 27r||/|p via Parseval's relation. With 
the outflow boundary condition (63) in rj space, it was shown by Eliasson (2001) for the onc-diniensional case that 
the energy integral 

||/||2= / r^' l/pdryrfx (66) 

is non-increasing in time, i.e. (i||/|p/(it < 0. In fact, due to that / is real valued, the symmetry condition /(x, —77, t) = 
f*{x,r],t) was used by Eliasson (2001), so that / is real-valued at 77 = 0. On the other hand, the total number of 
particles, 

N^ I 2TTfr,^ndx (67) 

is always conserved, and so is the total (kinetic plus electrostatic) energy of the system. 



W = / ( -TT 



drf 



^ \dx. (68) 

j)=0 ^ 



The decrease of ||/|p can be seen on as a loss of information about the finest details of the distribution function. 
Hence, the outflow boundary conditions in 77 space represents a dissipative process in which the highest oscillations 
in velocity space are removed from the system and a partial thermalization is allowed. In Fig. 5 we show the time 
evolution of the energy integral (66) relative to its initial value for a simulation of strongly nonlinear electrostatic 
waves, with the initial conditions (55) with amplitude A = 0.5 and wavenumber k^ = 0.5 (i.e. the same initial 
conditions as in Fig. 2), and using the outflow boundary conditions (63) at the boundary -q — rjmax = 30. Initially 
the energy integral decreases slowly, but at time i « 60 it decreases sharply. This time corresponds to "qmax/k, i.e. 
the time when the wave packet reaches the boundary 77 = rjmax as predicted by the solution of the reduced problem 
(62). The simulation was continued till t = 7000, and it could be observed that the energy integral decreased to a 
value of about 0.8130 after which it exhibited very small fluctuations (Eliasson, 2001). 

C. The relation between the outflow boundary condition and the Hilbert transform 

There is a simple relation between the outflow boundary condition (63) and the Hilbert transform, which in an 
infinite domain is defined as 

H[/](x) = -p r ^^dy (69) 

71" J -00 x-y 

where p denotes the Cauchy principal value. The outflow boundary condition (63) at the artificial boundary 77 = T^max, 
extended to an infinite spatial domain, is 

f^F-^H{k,)Ff, 77 = 77„ax (70) 

where the spatial forward and inverse spatial Fourier transforms are defined as 

00 

l/Ca; X 



F0= / 0(x)e-"==""=dx (71) 

and 

F-V - TT / '/'(fc)e''=="dfc, (72) 
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and the Heaviside function as 

The projection operator Q = F~^iJ(fcj.)F, acting on / as 

g[f]{x)=F-'H{k,)Ff (74) 

projects the function / onto the space of functions with only positive Fourier components in x space. In other words, 
the projection removes components with negative wavenumbers at the boundary and leaves components with positive 
wavenumbers unchanged. The Hilbcrt transform (69) has the property that 

H[e"===^](a;) = sign(fc,)ie"=^- (75) 

and it follows that the boundary operator Q can be expressed in terms of the Hilbert transform as 



6|/IW = i[/W-i«l/IW]=i 






(76) 



i.e., as an operator in real x space. Using Sohockij-Plemelj's formulas, we also have 



27r .5^0+ J_^ X — y + 10 

The integral formulations (76) or (77) of the boundary operator may open up the possibility to construct well-posed 
absorbing boundary conditions also for non-periodic problems, where the integrals are restricted to bounded limits. 

IV. THE NUMERICAL ASPECTS OF THE FOURIER TRANSFORMED VLASOV EQUATION 

We here discuss the numerical approximations used to solve the Fourier transformed Vlasov equation. As an 
example, we will discuss Fourier method for the one-dimensional Vlasov equation coupled with Gauss' law (Eliasson, 
2001) in some detail. Most of the results presented here carry over to the two- and three-dimensional cases (Eliasson, 
2002, 2003, 2007), and we therefore omit detailed discussions of the numerical methods in the multi-dimensional cases. 

A. Numerical approximations of the one-dimensional Vlasov equation 

We discretize the one-dimensional Vlasov equation on a rectangular, equidistant grid with the numerical box 
< a; < L in space, and < 77 < rymax in the Fourier transformed velocity space. The discrete function values at 
the grid points are enumerated such that f{xi, r]j,tk) ~ f^j with the spatial variable Xi = iAx, i = 0, 1, . . . , N^ — 1, 
and the Fourier transformed velocity variable rjj = jArj, j — 0, 1, ... , N^, where /S.x = L/N^, Aij — rj^a^/N^. The 
discrete time is obtained from the initial to = and then tk — t/t-i -I- Atk, k = \, 2, . . . , Nt The time step Atk may 
be kept fixed or varied dynamically to maintain numerical stability (Eliasson, 2002). 

The system (52)-(53) restricted to a finite domain, 

d? d^T 

-^ = i^^ - ir/^/, < 77 < r7„,ax, 0<x<L, (78) 

n 771 

— = l-27r/(a;,0,t), (79) 

are supplemented by the outflow boundary condition (63) at ry = 77max, 

/ - F-ii/(fc,)F/, r; = 77„ax, (80) 

and periodic boundary conditions in x space, 

/(x + L,r7,t) = 7(x,r7,t). (81) 
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Since / is real valued, we have the symmetry f{x,—rj,t) — f*{x,r],t), and hence the real part of / is even and 
the imaginary part odd with respect to rj (Armstrong et al., 1970; Eliasson, 2001). At 77 = one can therefore 
use a symmetry boundary condition, as discussed below. Equation (79) is solved numerically to obtain E, which is 
then inserted into right-hand side of Eq. (78); one can thus consider i<^ as a function of /. The outflow boundary 
condition at rj = rymax is performed by the boundary operator F^^H(kx)F, which removes all Fourier components 
with negative spatial wavenumbers (kx < 0) at the boundary. Discrete approximations arc used and for obtaining E 
in Eq. (79), for the boundary operator in Eq. (80), and for the rj and x derivatives the right-hand side of Eq. (78). 
The semi-discretized system can then be written as 



dfi.j 
dt 



Pifkj (82) 



where P is a grid function of all fij, representing the numerical approximation of the right-hand side of Eqs. (78). 
The unknowns fij are then advanced in time with the 4th-order Runge-Kutta algorithm: 

1. i^^V P(,/^), Vz,j- 

2. F^^f ^ P{f^ + p(i)At/2), yij 

3. F^y ^ P{j^ + F(2)At/2), \fij 

4. i^^y ^ P(,/^ + F(3)At), yi,j 

5. Jfr ^ft. + f (^? +2^^^ +2^? +^y). V^,J 
The steps needed for obtaining the approximation Pij are: 

1. Calculate the electric field numerically from Eq. (79). 

2. Apply numerically the boundary operator on the right-hand side, according to Eq. (80), for the points 
along the boundary 77 = 77max- 

3. Calculate a numerical approximation of the right-hand side of Eq. (78), for all points including the 
points along the boundary rj = rj^n^- 

Pseudo-spectral methods are used to approximate x derivatives and to integrate the Poisson equation, using the 
discrete Fourier transform (DFT) pair 

^^^W ^ '^^'"'^ '^''P ( "'^'''^iV" ) = DFT0(a;), (83) 

^t^i^j)^ J2 (t)ue^p[i2TrLu^^=DFT-^4>^. (84) 

uj=-{N^/2-l) 

The X derivatives are approximated in the pseudo-spectral method as 

^ « DFT-i [ik^BFTi^)] , (85) 

where kx — 2ttjuj/L, and the integration of the electric field in (79) is approximated by 

1 



DFT^i 



— DFT(1 - 2njf^„] 



(86) 



for all kx ^ 0, while the component of E corresponding to fcj. = is set equal to zero. 

In rj space, the derivative v = df/drj is calculated using the classical fourth order Padc scheme (Lele, 1992). For 
the inner points, the implicit approximation 

Vi,j-i + 4wij + Vij+i = — [fi,j+i - fi,j-ij , j == 1, 2, . . . , iV,, - 1 (87) 
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is used. At the boundary 77 = 0, the same approximation of the 77 derivative is used as for the inner points, taking 
into account the symmetry relations /i.-i = f*i and w^.-i = — Wj*i, 



or, for the real and imaginary parts. 



-<,i + 4wj,o + Wj,i = ^ ( fis - III 



(Ro) 



= 



2w 



(In 



,(Im) _ 



^Im) 



respectively. At the boundary 77 = r^max, the one-sided approximation 

is used. This gives a truncation error of order A?7'^ at the boundary. Equations (87), 
valued and one imaginary valued tri-diagonal equation system for each subscript i = 0, 1, 
Nri unknowns. 

At the boundary 77 = r/maxj the boundary condition (80) is applied with the approximation 

F-^H{k,)F<l){x,7j^,^,t) « DFT-i [i7(fc,)DFT((/)f „ 



(89) 
(90) 

(91) 

and (91) form one real 
Nx, each system having 



(92) 



where 



(X,Vn 



,t) is the right-hand side of (78) along the boundary rj = r^max and 0^jv its discrete approximation. 



B. The numerical representation of particle velocities 



3.0 



2.5 



VeffA?7 



Limit of exact effective velocity, v^s = v 



Maximum effective velocity, v^,™^ 




FIG. 6: The effective particle velocities produced by the 4th-order compact Fade difference scheme in 77 space. After Eliasson 
(2006). 

In order to investigate the impact of the difference scheme in the Fourier transformed velocity space on effective 
particle velocities, we study an approximation of the Fourier transformed Vlasov equation, in which the x and t 
derivatives are performed exactly, while the 77 derivative is approximated by Eq. (87). The difference scheme can 
formally be written as a difference operator Djj, giving rise to the difference-differential equation 



%-^liDj) + iE,f^O 



(93) 
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In order to study the effect of the difference approximation _D^ on the solution in v space, the definition (51) is inserted 
into Eq. (93), giving 



1 

2tt 



||^e""^-i|^2?,(e""')+ii?77/e""' 



du = (94) 



where the difference operator gives D^(e"'") — iWoff(w)e"'" with the effective particle velocity iwcff(w), and the term 
lErjfe^^'" gives rise to the term —e^^^^df /dv by an integration by parts, yielding 



1 

2^ 



dt ^ dx dv 



e'"" dv = 0, (95) 



and thus 



^+.e4-)^-i?^=0, (96) 

which is an approximation of the Vlasov equation (44), with the effective velocity V(,s{v) produced by the numerical 
differentiation in -q space. For the Fade scheme (87) we have 

3 sin(wA77) 

y^(v)^ ^ '/ (97) 

^ ' Ar] [2 + cos{vAr])] ^ ^ 

In the limit vArj — >■ 0, then Uoff — > v; see Fig. 6 where VesArj has been plotted as a function of vArj. The equations 
for the particle trajectories in {x,v) space, produced by Eq. (96), are 

. / N r / M 3 sin\v(t)Ar]] ,„„, 

x{t) - v.sHt)] = T^T^T— HtTTTTIT 98 

Arj {2 + cos[v(t)Ari\\ 

m^- — E[x{t),ti (99) 

Too 

where the dots in the left-hand sides denote time derivatives d/dt. Thus the particles are transported in x space with 
the effective velocity Vcs, which is periodic in v. If the product vArj is too large, then the approximation Vee ~ v 
breaks down; see Fig. 6. A maximum effective velocity, v^g^^ Arj = V^ ~ 1.73, can be found for vArj = 27r/3 « 2.09. 
It means that even though the largest represented velocity is given by the Nyquist limit Umax = tt/Atj w 3.14/Ar7, the 
highest effective velocity for transport of particles in x space is v^'^^^ — \/3/Arj. In numerical experiments, one has 
to choose a small enough Ary, so that important phenomena in velocity space, for example beams of particles, are well 
resolved, i.e., the velocities of these particles must fulfil v < v^^^^ with some margin. In the simulation performed to 
produce the results in the present section, particles were accelerated to velocities somewhat less than w = 5; see the 
left-hand panel of Fig. 3. The grid size was Arj = (120/600) = 0.2, which gives that vArj « 1.0 for these particles. 
According to the diagram in Fig. 6, the effective velocity is close to the limit of exact velocity for the value vArj = 1.0, 
and thus the particle velocities for the fastest particles are well resolved. The maximum effective velocity produced 
in the simulation was v^'^ — VS/Arj sa 8.6. 

C. The choice of domain and grid sizes 

When using the numerical algorithm to solve physical problems, it is important to know what is the computational 
domain and resolution in the real velocity space, i.e., what is the maximum velocity component fmax, used in the real 
velocity space to resolve the particle distribution function and what is the grid size Av used in the numerical solution. 
For given 77max and A77, the maximum represented velocity and the grid size in velocity space are given by 



and 



a;, '""" 



Av^ , (101) 

^max 



respectively. In view of the results in Fig. 6, one should choose A77 small enough such that the maximum velocity 
component Vmax is more than twice the maximum particle velocity one wants to resolve in the numerical solution, in 
order to avoid dispersive errors on the particle velocities. One must also choose ?7niax large enough so that fine enough 
structures in the velocity distribution function is resolved. 
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D. One-dimensional simulations of electron and ion holes 



In simulations of processes with several timescales, it is important that the numerical scheme does not introduce 
artificial heating of the electrons. Such effects can be a problem if a diffusion operator is introduced in velocity space 
to minimize the filamentation effects in the Vlasov equation. In the Fourier method used here, there is minimal 
heating effects of the electrons since only the highest Fourier modes in velocity space are absorbed by the outflow 
boundary condition in the Fourier transformed velocity space. As examples we will study the dynamics of electron 
and ion holes in an electron-ion plasma, which is governed by the Vlasov-Poisson system of equations for electrons 
and ions, 



dt 

dfi , df, 

■ V 



dx nip dx dv 



dt dx 

92$ _ e 

dx"^ Eg 



d^df. 



0, 



= 0, 



m,; dx dv 

{Ii{x,v,t) 



fe{x,V,t))dv 



(102) 
(103) 
(104) 



and where both the electron and ion dynamics are important. As initial conditions for the simulations, we have used 
Schamel's quasi-stationary solutions for electron and ion holes (Bujarbarua and Schamel, 1981; Schamel, 1971, 1979, 
1986). 
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FIG. 7: The electron distribution function for two electron holes at time t = (upper left panel), t = 155 ujp^ (upper right 
panel), t — 175ujp^ (middle left panel), t = 25lLOp^ (middle right panel), t = AGlujp^ (lower left panel) and t — 576 (lower 
right panel). The color bars go from dark blue (small values) to dark red (large values). After Eliasson and Shukla (2004a). 



Figures 8 and 7 show a simulation of interacting electron holes in and electron-ion plasma with mobile ions and 
realistic ion to electron mass ratio of 29500 for oxygen ions (Eliasson and Shukla, 2004a). Two large-amplitude 
electron holes were initially placed at x = — 40r£)e and x — AOrue, displayed at t = in the phase space number 
density plots in Fig. 7 (top left panel). A small electron density perturbation, centered at the two electron holes, 
were introduced as a seed for any instability. The ion density was initially taken to be homogeneous. We see in Fig. 
7 that at t = 155a;~g^ (top right panel), the two electron holes having started moving towards each other, and that at 
t = 175 ujpp they are colliding and are merging to a new electron hole (middle left panel), and that the newly created 
electron hole a.t t = 461 0;"^^ has propagated to x = SOroe (middle right panel), and at a later stage it has propagated 
to x = — 30r£)e where it remains throughout the simulation (lower panels). The reason for this complicated behavior 
of the electron holes can be understood by studying the interaction with the ions in detail. Due to their positive 
potential, the electron holes expel the ions and create local ion density cavities, which in turn eject and accelerate 
the electron holes. This can clearly be seen in Fig. 8, where the two electron holes start propagating at t •■ 



100 w-i 



and t Ki 130 Wpg, respectively. At t « 170 w^g, the two electron holes collide and merge into a new electron hole with 
larger amplitude, which propagates slightly in the positive x direction, and becomes trapped at a local ion density 
maximum at x sa 30r£ie; see the upper right panel of Fig. 8 for the ion density and the lower right panel for the 
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FIG. 8: The time and space evolution of the electron density (upper left panel), the ion density (upper right panel), the electric 
field (lower left panel) and the potential (lower right panel) associated with the two electron holes in Fig. 7. After Eliasson and 
Shukla (2004a). 



potential. After t ~ 400 a; g^, a new ion density cavity is created where the electron hole is centered, and at this time 
the electron hole is again ejected and accelerated in the negative x direction. At t w 480 w"^, the moving electron 
hole again encounters an ion density maximum located at x w — SOrue, where it is again trapped, performing large 
oscillations. We see that the electron holes remain stable during the acceleration by ion density cavities and survive 
on an ion time scale, much longer than the electron plasma period. 

Figure 9 displays the features of the ion and electron distribution functions for two colliding ion holes, where initially 
(upper panels) the left ion hole propagates with the speed uq = 0.9 VTi where vti — ykBTi/rrii is the ion thermal 
speed, and the right ion hole is standing. The ion and electron distribution functions associated with the ion holes are 
shown before collision at times t = Ocj"^ (upper panels) and t — 35.9 oj"^ (middle panels), and after collision at time 

t = 133 Wpj (lower panels), where ojpi — (noe^/eomi)^'^ is the ion plasma frequency. Figure 9 exhibits that the ion 
holes undergo collisions without being destroyed; thus they are robust structures. As can be seen in the right panels of 
Fig 9, the electrons have a non-Maxwellian, flat topped distribution in the region between the ion holes after collision 
has taken place. The velocity distribution function is plotted as a function of v/vtc in Fig. 10 at a; = 8.0 r^e- We 
see that the initial Maxwellian distribution (the upper panel) changes to a distribution with beams at u w ±0.6 vte 
(the middle panel) slightly before collision, and to a flat-top distribution with two maxima after collision (the lower 
panel). The reason for this phenomenon is that the two ion holes are associated with negative electrostatic potentials, 
and the electrons entering the region between the ion holes after collision must have a large enough kinetic energy 
to cross the potential barriers that are set up by the ion holes. Therefore, the region between the ion holes becomes 
excavated of low energy electrons. Also here the features of the electron distribution function survives on the ion 
time-scale, much longer than the electron plasma period. 



V. THE TWO-DIMENSIONAL VLASOV EQUATION 



We here discuss the extension of the Fourier transform method to the Fourier transformed Vlasov equation in two 
spatial and two velocity dimensions, including external and self-consistent magnetic fields. In the design of well- 
posed absorbing boundary conditions in the Fourier transformed velocity space, special care has to be taken with the 
magnetic field, which enters into the formulation of the boundary condition. 
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FIG. 9: The distribution function for the ions (left panels) and electrons (right panels) of two colliding ion holes, before the 
collision at times t = OtOpi (upper panel) and t = 35.9 LUpi (middle panel), and after the collision at f = 133cjpi (lower panel). 
The color bars go from dark blue (small values) to dark red (large values). After Eliasson and Shukla (2004b). 



Electron velocity distribution 




FIG. 10: The electron velocity distribution at x = 8ro, for t = Ouj ■ (upper panel), t = 35.9 lo^ (middle panel) and 
t = 133 tjp/ (lower panel). After Eliasson and Shukla (2004b). 



A. The two-dimensional Vlasov-Poisson system 
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The two-dimensional Vlasov-Poisson system for electrons reads 



df ^ df ^ df 
dt dxi 8x2 


e / 


axi o'a;2 


/92$ 92$X g 


no - 






-^2 t; h BV2 T, Bvi -— 

OV2 OVi OV2 



= 



CO /-OO 



/ dT7i d7;2 



CO ■/ —CO 



(105) 
(106) 
(107) 



where uq is the neutralizing heavy ion density background, fixed uniformly in space. The external magnetic field 
B(a;i,X2,t) is here directed along the x^, axis, perpendicular to the (a;i,X2) plane, and the electrostatic potential is 
calculated self-consistently from Poisson's equation. 
By using the Fourier transform pair 



f{xi,X2,Vi,V2,t) ^ 
f{xi,X2,m,V2,t) = 



CO /'CO 



/(:ri, X2, r;i, 772, i)e-'(''^''^+''^''^) dr72 dr,. 



— CO J —CO 

1 



(27r)2 



CO /'CO 



/(xi, X2, z^i, V2, i)e'(''i''i+''^"^) dz;2 d?;! 



CO •/ —CO 



the system (105)-(107) is transformed into 
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The systems (105)-(107) and (109)-(111) can be cast into dimensionless form 
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respectively, where time t has been normalized by uj I the velocity variables f 1 and W2 by fth,e: the Fourier transformed 
velocity variables m ^^d m by fj^,,' the spatial variables xi and X2 by tdo, the Fourier transformed distribution 
function / by no, the function / by JT^o^thc' ^he electric field components Ei and E2 by ^^th o^'oc'^e/e, the electric 
potential $ by v^^^ ^(m/e), and the magnetic field B by uj^^mc/e. 

In order to adapt the system (115)-(117) for numerical simulations, the computational domain is restricted 
to < a;i < Li, < a;2 < L2, < m < '?max,i and -r7i„ax,2 < m < ''7max,2- For negative 771, the symmetry 
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f{xi,X2,'rii,r]2,t) — f*{xi,X2,—'rii,—r]2,t) is used to obtain function values, if needed, owing to that the original 
distribution function /(x, v, t) is real-valued. It is therefore not necessary to represent the solution for negative 771 on 
the numerical grid. In the Xi and X2 directions, the periodic boundary conditions 

f{xi + Li,X2,'qi,V2,t) = f{xi,X2,m,'il2,t) (118) 

and 

/(xi,X2 + L2,rii^V2,t) = /(xi,X2,?7i,772,i) (119) 

respectively, are used. 

B. Outflow boundary conditions in Fourier transformed velocity space 

Similar to the one-dimensional Vlasov-Poisson system discussed above, we wish to design absorbing artificial bound- 
ary conditions in the Fourier transformed velocity space, so that the highest oscillations in velocity space can be cap- 
tured at the boundary and removed from the calculation. Hence, at the boundaries at 771 = ?7max,i and 772 = ±77max,2 7 
the strategy is to let outgoing waves pass over the boundaries, and to set incoming waves to zero. We explore the 
idea by studying the reduced initial value problem with only a constant magnetic field B = Bq, 

S - i^r^ - 'jr4r + ^o'^i ^ - ^o^^ ^ - o (120) 

ot oxiorji 0X20ri2 orj2 orji 

/(a;i,a:;2, 771,772,0) = /o(a;i,X2, ?7i, 772) (121) 

A Fourier transform in space {d/dxi — > ikxi and d/dx2 — > ikx2) gives a new differential equation for the unknown 
function f{kxi,kx2,m,'n2,t), 

df df df 

-^ + {kxi - Bo??2)^ + (fc.2 + Bom)^-^ = (122) 

I{kxi,kx2,mi'n2,t)t=o = foikxi,kx2,Vi,V2) (123) 

This is a hyperbolic equation for which the initial values arc transported along the characteristic curves, given by 



dt 

dV2{t) 
dt 



kxi-Bom{t) (124) 

kx2-Bomit) (125) 



Along the boundary ryi = 77max.ij Eq- (124) describes an outflow of data when kxi — -Bo?72 > and an inflow of data 
when kxi — So ^2 < 0. A well-posed boundary condition is to set the infiow to zero at the boundary, i.e., 

/m=r,_x,l - 0, kxi - B0V2 < (126) 

which can be expressed with the help of the Heaviside step function H as 

/ = H{kxl - -80772)/, 771 = ?7max,l (127) 

where 

The boundary condition (127) allows outgoing waves to pass over the boundary and to be removed, while incoming 
waves are set to zero; the removal of the outgoing waves corresponds to the losing of information about the finest 
structures in velocity space. 

Inverse Fourier transforming Eq. (128) then gives the boundary condition for the original problem (120) as 

f^F^^H{kxi-Bom)Fif, m=^max,l (129) 
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The operator Fj^ H{kxi — Bori2)Fi is a projection operator which removes incoming waves at the boundary 771 = ?7max,i- 
Similarly, the boundary conditions at 772 = ±77niax,2 becomes 



f^F^'H{kx2 + Bofii)F2f, 



m = ?7max,2 



(130) 



and 



/ = F^iH(-fc,2-Sor7i)F2/, 



m 



"'7max,2 



(131) 



respectively. 

In order to find well-posed boundary conditions in the 771 and 772 directions in the case when B = B{xi, X2, t) varies 
in time t and in the xi and X2 space with the periodicities Li and ^2, respectively, Eq. (115) is rewritten in an 
equivalent form as 



^-i^iA 
dt drji 



^-i.2i^or)(/V)-i^2^ 
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where we introduced the spatially averaged magnetic fields 



Boi{x2,t) = 
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Bdxi, 



Bo2{xi,t) = 
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Bdx2 



and the phase factors 



ii = exp 



irj2 / (B--Boi)dxi 
"'0 



h = exp 



rX2 

m / {B - B02) da;2 
Jo 



(132) 



(133) 



(134) 



We note that if B is periodic and continuous in the xi and X2 directions, then the integrals L'^{B — BQi)dxi and 
Jq^{B — i?02)dx2 are also periodic and continuous functions, and the xi and X2 derivatives can be approximated 
accurately by using the pseudo-spectral method. 

By studying the flow of data in the 771 direction for the function 0i = f9^ and the flow of data in the 772 direction 
for the function (j)2 = fO^ , one can find outflow boundary condition in the rji and 772 directions, similar to the 
conditions (129)-(131), as 
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/ 
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0iFr'ii(fcxi-772Boi)Fi(./0r'), 

92F^'H{kx2 + ViBo2)F2{fe2'), 
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m ^ ??max,l 
m ^ ?7max,2 
V2 = -J?max,2 



(135) 

(136) 
(137) 



respectively. In the case when B is independent of xi and X2, the boundary conditions (135)-(137) reduce to the 
conditions (129)-(131). The boundary operators are projection operators, which allow outgoing waves to pass over 
the T] boundaries of the domain and to be removed from the domain, while incoming waves are set to zero. The well- 
posedness of these boundary conditions was proven by showing that a positive deflnite energy integral is non-increasing 
in time (Eliasson, 2002). 



C. Electron Bernstein and upper hybrid Avaves in magnetized plasmas 



We will here give examples of some effects related to linear electron oscillations in magnetized plasmas. These 
include electron Bernstein modes that are exactly undamped according to Landau theory. Due to this fact, there 
is a recurrence effect in a weakly magnetized plasma, namely, there appears that waves can be periodically Landau 
damped only to recur in the plasma at a later time. We here compare simulation results with linear theory, in the 
form of dispersion corves for Bernstein mode waves and time-dependent analytic solutions of the Vlasov equation for 
a magnetized electron-ion plasma with immobile ions. 

The dispersion relation for the linear upper hybrid and electron Bernstein modes in a Maxwellian plasma is given 
by Crawford and Tataronis (1965) as 
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^pc 



exp(— A) 
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(138) 
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FIG. 11: Dispersion diagram and power spectrum (decibel) for electrostatic electron Bernstein (EB) and upper hybrid (UH) 
waves; iJuh = 4a;ce- After Eliasson (2008). 
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(139) 
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(141) 

Solving for oj' in Eq. (138) gives the relation between lo and k. Figure 11(a) shows the dispersion curves for the case 
'^uh — 4wco, and Fig. 11(b) shows the power spectrum in space and time from a simulation of the Vlasov equation 
with the same parameters. In the simulation, the obtained time series of the electric field component Ei was Fourier 
transformed in xi space and in time (using a Hamming window), to produce a power spectrum in a logarithmic scale. 
As can be seen in Fig. 11(b), the electrostatic energy is concentrated at the linear Bernstein eigen- modes, in good 
agreement with theory. In the long wavelength limit kro ^ 1, Eq. (138) reduces to the dispersion relation 



^uh 



3?)2 1,2 



(142) 



where Wuh = ('^pe+'-^co)^ i^ the upper hybrid frequency. However, taking electromagnetic effects into account, there 
are corrections in the long wavelength limit where the upper hybrid waves go over to the electromagnetic Z mode 
waves which we will discuss below. A zero-frequency (w = 0) mode, which is not a solution of the dispersion relation 
(138), can be seen in the power spectrum in Fig. 11(b); this "convective mode" has earlier been observed in numerical 
PIC simulations by Kamimura et al. (1978), and theoretically by Sukhorukov and Stubbe (1997). In terms of Landau 
theory, this mode is related to a pole in the initial condition, and not to a solution of the dispersion relation (138). 

The simulation was restricted to one spatial dimension, along the xi axis, where the simulation domain was 
set to < xi/r^, < AOn, < rjiVt^c < 15 and < ?72Wth.c < 15 and the number of intervals N^i = 50, N,^^ — 60 and 



2N„ 



120, respectively. The initial condition was set to 



f{xi,X2,m,V2,0) =n{xi)fo{T]i,r]2) 
where the perturbed relative density was set to a sum of waves with all possible wavenumbers. 



n{xi) = 



24 



1 + A'^ii sin(0.05ii2;i/rD) 



(143) 



(144) 



with the amplitude set to A = 0.0001, giving an almost linear problem. The Fourier transformed velocity distribution 
was set to a Maxwellian, 



h{vi,m) = 



no 

{2ny 
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(145) 
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In velocity v space, the Maxwellian would be 



no 1 



foivi,V2) = ::^— T^z^xp 



<,e Stt 






{vl + vD 



(146) 



The external magnetic field was kept constant in the simulation, with the ratio ujpc/uJcc = V 15 (giving uj^h = 4a;, 
The number of time steps taken was Nt = 42 530 and the end time iond = l^STWpJ^ with a fixed timcstcp. 
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FIG. 12: Pseudo-periodically damping and recurrence of electrostatic waves propagating perpendicularly to the magnetic field 
in a weakly magnetized plasma; the so-called Bernstein-Landau paradox, (a): Simulation results for different magnetic fields, 
(b): Analytic solutions by Sukhorukov and Stubbc (1997). 

According to linear theory, wave solutions of the unmagnetized Vlasov equation exhibit collision-less damping, 
while in magnetized plasma waves propagating perpendicularly to the magnetic field Boxt are exactly undamped, no 
matter how weak the magnetic field is. This is the so-called Bernstein-Landau paradox; it seems that the theory 
for magnetized plasma does not converge smoothly to the theory for unmagnetized plasma when the magnetic field 
strength decreases. 

This problem was investigated theoretically by Sukhorukov and Stubbe (1997) who derived analytical solutions of 
the problem. In Fig. 12, we have compared a numerical solution of the Vlasov equation with the analytic solution 
of Sukhorukov and Stubbe (1997) for a wave with the wavenumber k^ = 0.4 r^ and with different values on the 
magnetic field, such that ujcc/'-^pc =^ 0.4/10, 0.4/7 and 0.4/4, where Wco = e^oxt/"^c- The numerical result obtained 
from the Vlasov simulation can be seen in Fig. 12(a) show excellent agreement with the analytic results of Sukhorukov 
and Stubbe (1997) in Fig. 12(b). For the numerical results, where the electric field (normalized by its initial value) at 
X = is plotted as a function of time. The upper panels show the case with the weakest magnetic field -Boxt and the 
lower panel the strongest magnetic field. The horizontal time axis is scaled so that the tick marks are placed at each 
gyro period; in the upper panel the gyro period is igyro = (25 x 2Tr)u!~^, in the middle panel tgyro — (17.5 x 2Tr)u!~^ 
and in the lower panel tgyro = (10 x 2Tr)ujZj^. As can be seen in Fig. 12, the waves exhibit damping within the first 
gyro period, followed by a recurrence of the wave, which is again damped, etc, in an increasingly irregular pattern. 
In the upper panels of Fig. 12, with the weakest magnetic field, the field has the time to perform the largest number 
of oscillations within each gyro period and the electric field is also strongest damped before recurring at each gyro 
period. The "paradox" is resolved in the following manner: The waves exhibit damping within the first gyro period, 
given by tgyro = 27r/a;ce, where Wcc = ei?ext/we, and then the wave recurs the first time, followed by a new damping, 
et.c. In the limit of a vanishing magnetic field, -Boxt ~^ Oj it follows that cjco -^ and the gyro-period goes to infinity. 



U 



oo, and hence the wave will be damped, since it will take an infinite time for the first recurrence to occur. 



In the present numerical simulation, the initial condition was set to 

f{xi,X2,m,V2,0) =n(xi)fo{m,V2) 
with the relative perturbation 

n{xi) = [1 + As[n{kxiXi)] 
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(147) 



(148) 



and with the x-component of the wave vector set to k^i ~ 0.4 r^^. The amplitude of the wave was set to ^ = 0.0001, 
making the problem close to linear. The Fourier transformed velocity distribution was set to be a Maxwellian function. 



f(Tli-,m) == 



no 



exp 



- 2'^vl + vi)vth.c 



(149) 



The domain was set to < xi/ro < 27r/0.4, < ?7iWtho l£ 30 and —30 < ?72Wtho ^ 30, and the number of intervals 
N^i = 4, N,,, = Nri, = 150. 



D. Electromagnetic waves perpendicular to the magnetic field lines 



We restrict the Vlasov-Maxwell system to two spatial and two velocity dimensions, where the particles move in 
the (a;i,a:2) plane, the electric field E(a;i,a;2,t) is directed in the (a;i,a;2) plane and the magnetic fields B^xt £^nd 
B(a;i,X2,i) are directed in the 2:3 direction, perpendicular to the (xi,a;2) plane. We use the same normalization of 
variables as (Eliasson, 2003), i.e., t is normalized by w"^, vi and V2 by wth.o, xi and X2 by td, 771 and 772 by w^g, /c 
and /i by no, /o and /i by n^v^^^^, Ei and E2 by w^\ gro^We/e), $ by v^^^^{me/e), Soxt and B by Wpc(mo/e), Ai 
and A2 by Uth.o'Tic/e, Fi and r2 by ujpgVth,enie/^, P by nge, ji and j2 by Wth,e'T-oe, in which the Fourier transformed 
Vlasov equation for the ions and electrons attain the dimensionless form 
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respectively. The electromagnetic wave equations (27)-(28) take the form 
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and the first-order system (31)-(32) takes the form 
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The electric and magnetic fields are calculated as 
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and 
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respectively. The charge and current densities are calculated from the ion and electron distribution functions as 



p = (27r)2(/i - /e),,=,,=o - {27:f mfO ^ 3*(/c)),,=,.=o 
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respectively, and where S and 5ft denote the imaginary and real parts parts, respectively; the last equalities in 
Eqs. (159)-(161) follow from the symmetry condition where the real parts of the distribution functions are even and 
the imaginary parts are odd with respect to t]. 

We here consider electromagnetic waves that are propagating perpendicularly to the magnetic field lines. The 
waves have the electric field component in the plane perpendicular to the magnetic field, while the wave magnetic 
field is parallel to the external magnetic field. This configuration support the high-frequency X mode waves (X stands 
for "extraordinary") that can also propagate in vacuum as light, and the Z mode waves that connect to the upper 
hybrid resonance at short wavelengths shorter than the electron inertia length Ae — c/cope- In the low- frequency 
regime we have magnetosonic waves that connect to the lower hybrid resonance at wavelengths shorter than the 
electron inertial length. The O mode wave (O stands for "ordinary") has a component of the electric field along the 
background magnetic field direction, and cannot be simulated in the two-dimensional model discussed here, since it 
needs electron dynamics also along the magnetic field lines. The theoretical predictions for high- and low-frequency 
waves are compared with a Vlasov simulation in one spatial dimension (along the xi direction) and two velocity 
dimensions, of waves propagating perpendicularly to the external magnetic field. 

In Fig. 13(a) we show a comparison between dispersion curves obtained from the dispersion relation of cold fluid 
theory (Goldston and Rutherford, 1997) 



w2 u;2(w2-u;2|J 

and a Vlasov simulation, for the case Wuh/wco — 4 from which it follows that Wpo/wcc ■ 

used the ratio c/fth,e = 50 between the speed of light and the electron thermal speed. For large k, the we see that 

the fast X mode approaches the speed of light, while the Z mode wave approaches the upper hybrid oscillation, with 



(162) 



15. In the simulation, we 



frequency lj 



'^uh 



for large k. In the short wave length limit (very large k), thermal and kinetic effects 



are important, and the Z mode wave goes smoothly over to one of the upper hybrid and one of the electron Bernstein 
waves. The energy for the high frequency waves in Fig. 13(b) is concentrated at the linear dispersion curves for the 
fast and slow X modes, displayed in Fig. 13(a), in good agreement with theory. In Fig. 13(b), one can also see some 
weakly excited waves at the gyro harmonics uj/uJcc = 1, 2, 3, 4, which are waves not covered by the dispersion curves 
in Fig. 13(a), obtained from the cold plasma fluid model. The weakly excited oj/uJcc = 1 mode is an electromagnetic 
effect (Puri et ai, 1973) which can not be seen in the electrostatic case shown in Fig. 11 on page 24. 

In Fig. 14, we compare theory with a closeup of the low-frequency part of the energy spectrum obtained in the 
Vlasov simulation. For the low frequency electromagnetic waves perpendicular to the magnetic field. An approximate 
dispersion relation obtained from a cold fiuid description of ions and electrons is given by Goldston and Rutherford 
(1997) as 
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UJ^-UJ 



iH 



(163) 



with c/fth.o = 50, cjuh/'^cc = 4 and ujpc/uJcc = v 15. For numerical efficiency we use the mass ratio m^/mf, = 400 
between the ion and electron masses, which gives the ratios Wpi/wco = vT5/20 and ujd/ojcc = 1/400. Eq. (163) is 
solved for co/coce and displayed in Fig. 14(a). For large kxi, the dispersion curve approaches the lower hybrid frequency 
wih, approximately given by 



-"ih 



'^pi^ + (^ci^cc) ^ 



(164) 



and which is indicated in Fig. 14(a). For very small kxiXe < 1, we see in Fig. 14(a) that the dispersion curve 
approaches the one for Alfven waves, governed by the dispersion relation 
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FIG. 13: Dispersion diagram for the high frequency electromagnetic extraordinary mode, obtained from cold plasma fluid theory, 
and power spectrum (decibel) of the transverse part E2 of the electric field obtained from Vlasov simulation; cjuh = 4a;ce- After 
Ehasson (2003). 
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FIG. 14: Dispersion diagram for the ion extraordinary wave, obtained from cold plasma fluid theory, and power spectrum 
(decibel) of the transverse part E2 of the electric field obtained from the Vlasov simulation. The waves go over to Alfven waves 
at small k^i and to lower hybrid hybrid oscillations at large kxi; uJu\i = 4tJce, mi/rrie = 400. After Eliasson (2003). 



where va = o-^d/'-^pi is the Alfven speed. The energy spectrum for the low frequency waves in Fig. 14(b) shows good 
similarity with the dispersion curve for the low frequency wave in Fig. 14(a). The width of the energy bands in the 
power spectrum is the frequency resolution obtained in the simulation; a longer simulation would resolve the waves 
more. The frequencies of the waves in Fig. 14(b) for large k is slightly higher than the corresponding frequencies in 
the dispersion diagram in Fig. 14(a), which probably is a thermal effect, not included in the cold plasma fluid model. 
In the Vlasov simulation presented here, the simulation domain was restricted to one spatial dimension with 
< xi < 20007r and two velocity dimensions, plus time. The Fourier transformed velocity domain was < 771 < 10 
and —10 < ?72 ^ 10 for the electrons, and < 771 < 200 and —200 < 772 < 200 for the ions, with the number of intervals 
Nxi = 100 in space and iV^j^ = 30 and ^N^j^ = 60 in the Fourier transformed velocity space. The initial conditions 

for the electrons and ions were set to fc{xi,X2,r]i,r]2,0) = ?^(a;)/o,o(f7i, '72) and fi{xi,X2,r]ij''l2,0) = "(a;)/i,o(?7i,?72), 
respectively, with the density perturbation 
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n{x) = 1 + 10"^ Y^ h sin(0.05iia;i) 



(166) 
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so that all wavcmodes were excited with low amplitudes, making the problem close to linear. The electron and ion 
distribution functions were set to 
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respectively. The number of time steps taken was Nt — 97 280; the end time was tend — 8340. No numerical dissipation 
was used. The ion-electron mass ratio was mi/rric — 400 and the ion and electron temperatures were set equal, Tj = To, 
giving the factor (wth.e/wth.i)"^ = 1/400 in Eq. (168). 

VI. THE THREE-DIMENSIONAL VLASOV EQUATION 

We here discuss the extension of the Fourier transform method to the Fourier transformed Vlasov equation in full 
three spatial and three velocity dimensions, including external and self-consistent magnetic fields. Again, in the design 
of well-posed absorbing boundary conditions in the Fourier transformed velocity space, special care has to be taken 
with the magnetic field, which enters into the formulation of the boundary condition. 

The Fourier transformed Vlasov equation (13) can be cast into the normalized form 
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where t is normalized by uj^^, v by wxe, x by tdo, V by w^^\ fa by no, E by v^^r-^l{mc/e), and B by Wpc(mo/e). 
Here we defined Qa — qa/s and Ma — ma/rnc, where we assume electrons and singly charged ions, so that Qi — 1, 
Qe = ~lj M\ = "^i/'^oj a-iid Mo ~ 1. The equations for the potentials, (31) and (32), take normalized form 



dA_ 

dt ' 

and the electric and magnetic fields (33) and (34) arc obtained as 






-V^A+j, 



'^Te 



E 
B 



-Vx$-r, 

Bo + Vx X A, 



(170) 
(171) 



(172) 
(173) 



where $ is normalized by w^^(?7ic/e), A by WTo("T-c/e), F by u;pcWTc('™c/e). p by eno, and j by enowxe- Using Eqs. 
(14) and (15) in Eqs. (7) and (8), the normalized charge and current densities are obtained as 



and 



respectively. 



p=(27r)3[SR(/i)^=o-5ft(,/c)^=o], 



j = (27r)3[V^3(/i)^=o - V^3(/e)^=o], 



A. Restriction to a bounded domain 



(174) 



(175) 



In order to adapt the Fourier transformed Vlasov Maxwell system for numerical simulations, it must be restricted 
to a bounded domain. The computational domain is restricted to < xi < Li, < a;2 < i2j < 0:3 < Z/3, < ?7i < 
??max,ia, -'7max,2a < V2 < ??max,2a, and -T/max.Sa < V3 < ??max,3a- Hcrc a (cqual to 6 for elcctrons and i for ions) is 
introduced so that different domain sizes can be used in t] space for the ion and electron distribution functions. For 

negative rji, the symmetry f(xi,X2,X3, 771 , 772 , 773 , i) = f*{xi,X2, X3, — ?7i, —772, ^%ji) is used to obtain function values; 
it is therefore not necessary to numerically represent the solution for negative ryi if the solution is represented for 
negative 7^2 and 7^3. 
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B. Outflow boundary conditions in Fourier transformed velocity space 



In this section, wc will derive well-posed boundary conditions for the Vlasov equation in r] space. Writing out the 
terms of the Fourier transformed Vlasov equation (169), we have 



dt dxidrji dx2drj2 dx^drj^ Ma 



i{Eir]i + ^2*72 + E-im)fa 



dfa dfa 

+ {B2m - -83%)^ + {B3V1 - BiT]3)—^ + {Bir]2 - B21JI, 

drji dr]2 d-qz 

In position space, periodic boundary conditions 

fa{xi +Li,X2,X3,T]i,r]2,r]3,t) = fa{xi,X2,X3,r]i,r]2,V3,t), 
fa{xi,X2+L2,X3,T]i,r]2,m,t) = /a (^i , ^2 , X3 , ?7i , 772, % , t), 
fa{xi,X2,X3+L3,T]i,r]2,m,t) = /a (^i , ^2 , X3 , ?7i , 772, % , t). 



(176) 



= 0. 



(177) 
(178) 
(179) 



are used for both the distribution functions and the electromagnetic fields. The artificial boundaries at 771 = r;max,iaj 
772 = ±'7max,2Q and ?73 = ±77niax,3Q must be treated with care so that they do not give rise to reflections of waves or 
to instabilities. The strategy is to let outgoing waves pass over the boundaries, and to set incoming waves to zero. 
The problem of separating outgoing waves from incoming waves is solved by employing the spatial Fourier series 
expansions (transforms) . In order to explore the idea, one can study the reduced initial value problem with a constant 
external magnetic field B = Bg and a zero electric field E = 0, 



dt dxidrji dx2drj2 dx3dri3 M^ 



{B2,aV3 - B3^om 



dm 



(B3.0V1 - Bi^of]3)^ + {Bi^Qf]2 - B2.om)^ 

dr]2 dri3 



0, 



fa{xi,X2,X3,'qi,'n2,m,t)t=0 = /ao(a;i , 2:2 , X3 , ?7i , ?72, %) • 

By introducing the spatial Fourier series pairs in (xi, X2, X3) space, 

^1 Jo 



27rii 



9i,iie 



ik^-^ xi 



h,t2 = F202 = 
^2 = F-1^2 = 



2]^^ — 00 

1=0, ±1, ±2, 
L2 



L2 Jo 

00 

- E 



02(a;2)e-'^-^"Mx 



62,., 6"^-^^=^ 



27ri2 



, i2-0, ±1, ±2, 



(180) 
(181) 

(182) 
(183) 
(184) 
(185) 
(186) 
(187) 



and 



»3,»3 



^3'?^3 



L 



3 Jo 



03(a;3)e"'"^3"^dx3, 



i>3 = F3 (j)3 



E '^3,^3e''-^"^ 



kx3 = -,— , h = 0, ±1, ±2, 
^3 



(188) 
(189) 
(190) 
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respectively, and Fourier-transforming Eq. (180) in all spatial directions, one obtains a new differential equation for 
the unknown function /(fc^^, fc^^s^ ^2:37 ?7i, '?2, %, *), 



dt ^ 


f^Xi 




- B^^om) 


+ 


^X2 


-M?^-^- 


- Bi^om) 


+ 


f^X3 


" ji-{Bifiri2 - 

Ma 


- B2firii) 



dfa 

dm 

dfa 

dm 

dfa 

dm 



(191) 



0, 



with the initial condition 



faikxi , kx2 , kx3,m7V2, m, t)t=0 = faoikxi , ^xa , ^xg , J?!, J?2, %)• 



(192) 



Equation (191) is a hyperbolic equation for which the initial values are transported along the characteristic curves, 
given by 



-^r^^'^^-M-a^^'^''^'-^''''^'^' 

X2 j^{B3fiVi- Bi^oVs), 



dV2{t) 

dt 
dysjt) 

dt 









{Bifim - B2,om)- 



(193) 
(194) 
(195) 



Along the boundary rji = ?7max,Qii Eq. (193) describes an outflow of data when k^-^ — {Qa/ Ma){B2.om " ^3,0^2) > 
and an inflow of data when k^^ — (Qa/ Ma){B2.om ~ B^^om) < 0. A well-posed boundary condition is to set the 
inflow to zero at the boundary, i.e.. 



(ja)rii=ri„ 



Qo 



0, kx^ - j^{B2flm - B^fim) < 0: 



(196) 



which can be expressed with the help of the Heaviside step function H as 



U^H 



XI - J^iB2,0V3- B3^oV2) 



fa, m = Vn 



(197) 



where 



HiO = 



1, e>o 

0, ^ < 0, 



(198) 



The boundary condition (197) allows outgoing waves to pass over the boundary and to be removed, while incoming 
waves are set to zero; the removal of the outgoing waves corresponds to the losing of information about the finest 
structures in velocity space. Inverse Fourier transforming Eq. (197) then gives the boundary condition for the original 
problem (180) as 



U = Fl'H 



kxi - j^{B2fim - B^flm) 



Fl/a: m — '7max,al- 



(199) 



The operator Fj^ H[kx^ — {Qa/Ma)(B2.om ~ -B3.o?72)]Fi is a projection operator which removes incoming waves at 
the boundary m — '?max,ai- Similarly, the boundary conditions at m — ±?7max,a2 and 773 = ijymax.aa become 



fa -F2 H 



/a - F2 H 



kx2 - j^{B3,oVi - Bifim) 

Qc 



kx2 + j^{B3,om - Biflm 



F2/a, m — '7max,a2, 

F2./a, m — ^'7max,a2, 



(200) 
(201) 
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and 



/a - Fg H 



fa - F, H 



«a;3 - J^{Bl,0V2 ~ B2fim) 



kx3 + j^{Bi^ori2 - B2.Qr]i) 






(202) 
(203) 



respectively. 

In order to find well-posed boundary conditions in r) space in the case when B varies both in space and time, 
Eq. (176) is rewritten in the form 





■ . d 

dxi 


-Ai 


iUOai) 


dr]2 


■ . d 

3X2 


-I302 


{tO~2) 


dm 


■ . d 

dx3 


-I303 


iLo-J) 



(204) 



K-^l^l + ^2112 + E^mYfa = 0, 



where the phase factors O^i, 9a2 and 0^3 are 



^cti = exp 
9a2 = exp 



and 



respectively, and where 



ya3 = exp 



xi 

1 / {Pi~ /3oi)da;i 



1 / (^2-/3o2)da;2 





1 / (/33-/3o3)da;3 





Qo 
Qo 



(-B2?73 - -83^2) 



and 



p2 = j^iBa-qi- BiTjs) 



Pa = Jf(^^^^ - B2V1), 



/9oi 

/302 



Li 



L 



Pi Axi 



f32dx2 



2 Jo 



/3o3 = ^— / Psdxs. 

^3 Jo 



(205) 
(206) 

(207) 

(208) 
(209) 

(210) 



The form (204) of the Vlasov equation makes it possible to introduce stable numerical boundary conditions in t/ space 
in a systematic manner. Furthermore, Jq^ {Pi — /3oi)dxi, /g^^(/32 — /3o2)dx2 and Jq^IPs — Poa)'^^^ are continuous and 
periodic in x space if B is continuous and periodic in x; this is the reason for the subtraction of the mean values /3oi, 
/3o2 and /3o3 in the integrals. 

By studying the flow of data in the 771, 772 and m directions for fa&ai j /a^a2 ^'^d fada3 > respectively, one finds the 
outflow boundary conditions to be 



fa — Sal^i H{kxi ~ Poi)^l{fa()al )' Vl = Vmax,al, 

fa — ^a2F2 H{kx2 ~ /3o2)F2 (/a6'^2 ) ' ^2 = '7max,a2, 

fa — ^a2F2 H{ — kx2 + /3o2)F2 (/Q^ct2 ) ' V2 = — ?7max,Q2, 

fa — dasFs H{kx3 — PosWsifada^), rj2 — ?7niax,a2, 

fa = OaiV3^H{-kx3 + /3o3)F3(/a^Q3 ), m = -r7max,a2- 



(211) 

(212) 
(213) 
(214) 
(215) 
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In the case when B is independent of x and t, the boundary conditions (211)-(212) reduce to the conditions (199)- 
(203). For the case where the domain is extended to negative rji, we also have the boundary condition 



L = 0^iF^'H{-k,, + (3oi)FiiUe-l), 



m 



^max.al ■) 



(216) 



which was used at a stage in the numerical algorithm. The well-posedness of these boundary conditions were proven 
by using that a positively definite energy integral was non-increasing (Eliasson, 2007). 



C. Electromagnetic electron waves 




FIG. 15: a) Amplitude spectrum (10-log scale) of E2 for waves propagating parallel to the external magnetic field (6 = 0). b) 
Dispersion relations for the high-frequency R and L waves and the low-frequency whistler wave, propagating along the external 
magnetic field {6 = 0). c) Amplitude spectrum of E2 for waves propagating with an angle 9 = 7r/4 to the external magnetic 
field, and d) dispersion relations for electron waves propagating with an angle of ^ = 7r/4 to the external magnetic field. Here 
oJco/t^Jpc = 1/4. After Eliasson (2007). 



The general dispersion relation for electron waves in a cold collisionless plasma with an external magnetic field is 
given by the Appleton-Hartree dispersion relation (Stix, 1992) 



a;2 



= 1- 



2(^2 






± Wr,oA 



(217) 



where A = [ 



^1^ sin^ I 



+ 4l^-2(^2 



'cos^^]^/^ and is the angle between the external magnetic field and the 



wave vector. In order to assess that the Vlasov code reproduces these well-known wave modes in the plasma, we have 
simulated electromagnetic waves propagating at different angles to the magnetic field; see Figs. 15 and 16. In the 
simulations, we restricted the problem to one spatial dimension, along the x\ axis, and used three velocity dimensions. 
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FIG. 16: a) Amplitude spectrum (10-log scale) ol Es for waves propagating perpendicularly to the external magnetic field 
{9 — 7r/2), with the electric field perpendicular to the external magnetic field, and b) dispersion relations for the high-frequency 
slow and fast extraordinary (X-) modes, propagating perpendicularly to the external magnetic field {9 — 7r/2). c) Amplitude 
spectrum of E2 for waves propagating perpendicularly to the external magnetic field {9 — 7r/2), with the electric field parallel 
to the magnetic field, and d) dispersion relations for the ordinary (0-) mode propagating perpendicularly {9 = it/2) to the 
magnetic field. Here Ucc/i^pc = 1/4. After Eliasson (2007). 



The initial condition for the electron distribution function was a Maxwellian distribution (in normalized units) 



/e = (27r) exp 



(218) 



while for the ions we used 



/i = (2^)-^cxp 



2 V Terrii 



1/2' 



(219) 



with Ti/Tc = 100 and nii/me = 1836, and we chose c/vtc = 50 for the electromagnetic waves. The magnetic field 
strength was chosen such that ujce/^pe = 1/4, i.e. in our scaled unit we have |Bo| = 1/4. A low-amplitude noise 
(random numbers) were added to the vector potential A and to T so that all wave modes in the system were excited. 
The numerical parameters were chosen as N^ 



40. 



dimensional units), 



Nx2 = ^2:3 ~ I, Li = n X 10'^ (corresponding to 207rc/a;po in 



^m = Nv2 



N, 



m 



20, 



^max,el — ^niax,c2 — ^max,c3 



= 6, and ?7max,il = ^max,i2 = ??i 



max, 13 



200. 



The simulations were run with 8000 time-steps with the fixed time interval At — 0.14. In Figs. 15 and 16, we have 
Fourier transformed the electric field in space and time (with a Gaussian time window) to obtain the spatio-temporal 
wave spectrum. In panel a) of Fig. 15, we show the power spectrum for the transversal electric field component E2, 
for waves propagating along the magnetic field lines. It is clearly seen that the wave energy is concentrated along 
the dispersion curves of the electromagnetic right-hand (R) and left-hand (L) circularly polarized waves, shown in in 
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panel b). They are given by the dispersion relation 



= 1 - 



1 - Wco/w 



(220) 



and 



,2j,2 



c^k 



po 



/-^ 



1 + Wco/w 



(221) 



respectively, obtained by setting = in Eq. (217). The R-wave is divided into a high-frequency branch (having the 
highest frequency) and the low-frequency electron whistler branch. We next made a simulation of waves propagating 
obliquely to the external magnetic field, which was chosen as (i?oij -B02, -B03) = (0.25, 0.25, 0)/-\/2. The resulting 
amplitude spectrum of E2 is presented in panel c) of Fig. 15, while the solutions of the dispersion relation (217) is 
plotted in panel d). Here, we can see the emergence of the slow X-mode which has a resonance somewhat higher 



than the plasma frequency uj 



UJr 



3. Comparing panel c) and d), we see that the wave energy is concentrated at the 



dispersion curves. In Fig. (16), we are considering waves propagating perpendicularly to the magnetic field. Here, 
the external magnetic field is given by (-Boi, ^02, ^03) — (0, 0.25, 0), and the energy spectrum in panel a) and c) are 
for the perpendicular (to the magnetic field direction) and parallel electric field components E^ and E2, respectively. 
The wave energy is concentrated at the dispersion curves for the cold plasma fast and slow X-modes displayed in 
panel b) and the 0-mode plotted in panel d). The cold plasma dispersion relation for the ordinary (0-) mode and 
extraordinary (X-) mode, perpendicular to the magnetic field, is given by 



and 



C2fc2 


= 1- 


W2 






W2 




C2fc2 


= 1- 


-^c 




(-^--^c) 


W2 


W2 


(c. 


2 - C.2 _ ^2J 



(222) 



(223) 



respectively, obtained by setting 6 =- -k jl in Eq. (217). Also seen in panel a) at (jj/oj-pe = 0.5 is an excitation of an 
electron Bernstein mode which starts at cj = 2L0ce at small wavenumbers and has a resonance at cj = lo^c for large 
wavenumbers. 



D. Temperature anisotropy driven wfhistler instability 



In magnetized plasmas, there are often different temperatures parallel and perpendicular to the magnetic field 
direction. In this case, we may have a firehose instability if Tgii > Te_L, or a whistler instability if Tc±_ > T^h. The 
latter case can have relevance both for the sun and for the Earth's magnetosheath (Gosling et at, 1989; Zhao et al., 
1996). In order to study the growth and saturation of the whistler instability, we carried out a simulation where 
the initial condition for the electrons was taken to be a bi-Maxwellian distribution function, where we used the 
temperature ratio Tg^/T^u = 8. In the Fourier transformed velocity variables used in the simulation, the electron and 
ion distribution function takes the form 



while for the ions we used 



/e = (27r)-3 



/i = (27r)-^cxp 



exp 



(772 + 8r,2+8r;2) 



ivl 



vl 



vl) 






1/2' 



(224) 



(225) 



with Ti/T(,|| = 100 and rai/nic = 1836. We use the same numerical parameters as in section VIC, except that 
'?max.c2 = 'ymax.cS — 3 and wc usc a higher resolution in space such that N^i = 80. An adaptive time step was used 
in the simulation to maintain numerical stability (Eliasson, 2007). The numerical results are displayed in Figs. 17 
and 18. The time-dependence of the maximum amplitude (over the spatial domain) of the perpendicular electric and 
magnetic field components E2 and B2 are shown in Fig. 17, and we see an exponential growth of the perpendicular 
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2000 



FIG. 17: a) The maximum amplitude of the paraUel and perpendicular electric field components E^ and E2, respectively, and 
b) the perpendicular magnetic field component B2 (The parallel electric field Ei is shown only for times t > 500 a;";?). After 
Ehasson (2007). 



electric field component £'2, with a growth rate of wi ~ O.OlTwpc [indicated by the solid line in panel a)] of both the 
electric and magnetic field. The initially almost purely electromagnetic waves saturate nonlinearly by exciting the 
electrostatic field component Ei (see the upper panel of Fig. 17), and the amplitude of the perpendicular magnetic 
field fiuctuations are at this point about 10 % of the external magnetic field. In order to compare the simulation result 
with theory, we have plotted the spatio-temporal amplitude spectrum of the perpendicular electric field component E2 
in panel a) of Fig. 18 (in a 10-log scale) where the Fourier transform in time was taken for waves between LUpct = and 

t = 600, with a Gaussian time window. In panel b), we have plotted the quantity wi — (l/ii)ln[|i?2(ii)|]+constant 
where E2 is the spatial Fourier transform of the electric field component E2 at time t = ti and ii = 700 w~^. Panel 
b) gives a rough estimate of the growth rate for different wavenumbers in the simulation. We see that there is a 
significant growth rate of waves with wavenumbers between fcc/wpe « 1.2 and kc/ujpc ~ 2.6. We have solved the 
dispersion relation for whistler waves in a plasma with a bi-Maxwellian electron distribution. In the one-dimensional 
case, and for immobile ions, it is (Stix, 1992) 



■^pe 



2^2 



fcnC- 



c_L 






- 1 



{OJ - Wcc)Tc_l/Tc|| +Wc. 

\/2A:||Wxc|| 



/2knVTo 



(226) 



where Wxey = {kBTe\\/me) ' is the parallel electron thermal speed and 



(227) 



is the plasma dispersion function. In the dispersion relation (226), we have neglected the electromagnetic displacement 
current term [corresponding to the first term in the right-hand side of Eq. (220)]. In panels c) and d) of Fig. 18, we 
have plotted the real and imaginary parts of the frequency, obtained from the dispersion relation (226), where we have 
used the simulation parameters Wce/wpo = 0.25, Tf.±/T^\\ = 8 and c/uTeii = 50. Comparing panels a) and c), we see 
that for the undamped waves at small wavenumbers, the wave energy of the waves in the simulation is located along 
the dispersion curve of the whistler wave. Panels b) and d) show that the spectrum of the growing waves matches 
approximately the waves with positive growth rate obtained from the dispersion relation (226). We note that the 



maximum growth rate oji 
Fig. 17. 



0.017 Wpe in panel d) of Fig. 18 agrees well with the measured growth rate in panel a) of 



37 



a) 11 

1 


• "IPW 


o.s 


E " "' 


o.s 






Q.7 






u D.6 






D.4 


? 


*. 


0,3 






Q.£ 


% 





n.1 






-5 


■4 -3 -Z -1 1 2 3 4 5 



c) 



I' 



03 














■ 


0^ 






^""----^^ 








■ 


0.1 
3 








\./ 






■ 


i 


-1] 


-2 





2 


4 


a 




FIG. 18: a) The spatio-temporal amplitude spectrum of the electric field component E2, and b) the spatial amplitude spectrum 
at time t = TOOcjpe. c) The real and d) imaginary parts of the frequency for whistler waves, obtained from the dispersion 
relation (226), in an anisotropic plasma with T±/T\\ = 8, c/vtc\\ = 50 and ujcc/^-pc ~ 0.25. After Eliasson (2007). 

VII. EXTENSIONS TO INCORPORATE RELATIVISTIC AND QUANTUM EFFECTS 

We here very briefly discuss the extensions of the Fourier technique for relativistic and quantum Vlasov equations. 

A. The relativistic Vlasov equation 
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FIG. 19: The function G{r]). 

In the relativistic Vlasov equation, the relativistic gamma factor comes into play for particles moving close to the 
speed of light. This complicates the use of the Fourier transform technique to solve the Vlasov equation. As a model 
example we consider the one-dimensional Vlasov-Poisson system 



df pdf df „ 
ot 7 ox op 



dE 
dx 



= 1- 



fdp 



(228) 
(229) 



where the relativistic gamma factor is 7 = ^Jl -t-p^, and / describes the distribution of electrons in {x,p) space 
where p is the momentum. Here the distribution function / has been normalized by no/c, time t by w"^, space x by 



Ae = cjijjpe-, niomcntum p by mgC, and the electric field E by me(? / eX^- Using the Fourier transform pair 

/oo 
^TT J-oo 

we have the Fourier transformed Vlasov-Poisson system 



f-*aS,(^*^")+^^^^"^0' 



dE 
dx 



1 -27r/(x,?7,t)^=o, 



where * denotes convolution over 77 space. Here 

1 roo 



G 



2tt 



TT 
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(230) 
(231) 

(232) 
(233) 

(234) 



where Kq is the Bessel function of second kind of order 0. The function G{rj), plotted in Fig. 19, grows like 
[— 7_E — ln(|??|/2)]/7r for jryj < 1, where ^e ~ 0.5772156649 is the Euler-Mascheroni constant, and falls off like 
exp(— |77|)(27r|?7|)~^'^ for \ri\ > 1, and has the property that /^ G drj = 1. The relativistic corrections are con- 
tained in G. For a weakly relativistic plasma, the distribution function / is much wider and smoother in 77 space 
than G, and hence we have G * / « /, i.e., G then has the property of Dirac's delta function and we retain the 
non-relativistic Fourier transformed Vlasov equation treated by Eliasson (2001). The numerical implementation of 
the convolution by G and the well-posed absorbing boundary conditions in 77 space are unsolved problems, and so are 
the extensions to higher dimensions. Since the function G in Fig. 19 falls of exponentially for large 77, the convolution 
integral in (232) can possibly be approximated by a truncated operator with compact support, and the problem with 
absorbing boundary conditions could potentially be solved using normal mode analysis similar to that of Engquist 
and Majda (1977, 1979). 

B. The Quantum Vlasov/Wigner equation 

The quantum analogue to the Vlasov-Poisson system is the Wigner-Poisson model (Markowich et al., 1990; Roos, 
1960; Tatarskii, 1983; Wigner, 1932). In three dimensions, the Wigner equation for electrons can be written 



df 



V V/ 



d^Xd^v'e^'^'^^-^'^^/^ 



dt ' ' ••' (27r)3n'* 

which is coupled with the Poisson equation with immobile ions 



x+^,t 



A 

X ,t 

2' 



/(x,v',t). 



VU = 



eo 



fdv — no 



One can show that the Wigner equation converges to the Vlasov equation in the formal limit ?i — > 0, 

at me c'v 



(235) 



(236) 



(237) 



It turns out that the Fourier technique in velocity space is well suited to solve the Wigner equation. As an example 
we will study the ID Wigner-Poisson system (Manfredi, 2002) 



dt dx 



iem„ 



(27r)r 



^irrie {v~v )X/h 



x+^,t 



x-^,t 



f{x,v',t)dXdv' 



(238) 



dx'^ 



eo 



fdv -no] . 



(239) 
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As it stands, the Wigner is awkward to solve numerically. However, introducing the Fourier transform pair in velocity 
space 



f{x,v,t) 



f{x,r,,t)e-'^^ dr,, 



(240) 



we obtain 



dt dxdrj h 



/(x,77,i) = — / f{x,v,t)e^^^ dv 



<i)\x + 



2m„ 



2m„ 



/-O, 



(241) 



(242) 



d) € 

-— = — [27r/(a;, ij, t)^=o - no] 



(243) 



which is simpler to solve numerically. By using a pseudo-spectral method in space, the spatial shifts by ±hri/2me in 
Eq. (242) is converted to mulitplications by exp{±ikfiri /2me) , where k is the wavenumber. In 77 space, we can apply 
the same absorbing boundary conditions. 



/==F-i[i/(fc)F/] at 77 = 77,, 



(244) 



as for the Vlasov equation. Hence the cxisitng Vlasov codes are easily modified to simulate the Wigner equation; see 
for example the work by Marklund et al. (2006) where the Wigner equation for broadband electromagnetic radiation 
in a plasma was solved with the Fourier method as described here. Finally we mention that other numerical methods 
for solving the Wigner equation exist in the litterature, for example operator splitting methods (Arnold and Ringhofcr, 
1996; Suh et al., 1991). 

VIII. CONCLUSIONS 

In this paper we have given a review of simulations of the Vlasov equation in higher dimensions. In this method, the 
Vlasov equation is Fourier transformed in velocity space, and the resulting equation is solved numerically. We have 
discussed the main difficulties solving the Vlasov equation with a grid-based solver, namely that in some problems, the 
particle distribution function becomes filamented in velocity space due to phase mixing. This can lead the recurrence 
of the initial condition on the numerical grid (the so-called recurrence phenomenon) , which in turn leads to unphysical 
oscillations and instabilities in the simulations. By designing outflow boundary conditions in the Fourier transformed 
velocity space, the highest Fourier modes in velocity space are allowed to propagate over the boundary and to be 
removed from the simulations. In that way an effective dissipation is allowed in the Vlasov equation, and the numerical 
recurrence phenomenon is strongly reduced. On the other hand, Fourier modes that have not reached the boundary 
in the Fourier transformed velocity space are not damped by the numerical method. In a sense, the method represents 
the minimum numerical viscosity possible one can introduce to the numerical simulation, which removes the recurrence 
phenomenon. The extension to multiple dimensions was also possible, by careful consideration of the magnetic field 
in the boundary conditions in the Fourier transformed velocity space. In this manner, the boundary conditions could 
be made strictly local in time and to only include boundary points in the Fourier transformed velocity space. The 
boundary conditions are highly absorbing and have been proved to be well-posed by energy estimates (Eliasson, 2001, 
2002, 2007). The method may be an attractive alternative to existing methods for solving the Wigner equation, which 
describes the evolution systems of quantum particles. For the unmagnctized Wigner equation for charged particles, 
the method is directly applicable with only minor modifications. Future developments of the Fourier method could 
include relativistic effects and the extension to magnetized quantum plasmas. 
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